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Abstract 

This paper presents a new joint inversion approacii to siiape-based inverse prob- 
lems. Given two sets of data from distinct physical models, the main objective is to 
obtain a unified characterization of inclusions within the spatial domain of the physi- 
cal properties to be reconstructed. Although our proposed method generally applies 
to many types of inversion problems, the main motivation here is to characterize 
subsurface contaminant source-zones by processing down gradient hydrological data 
and cross-gradient electrical resistance tomography (ERT) observations. Inspired 
by Newtons method for multi-objective optimization, we present an iterative inver- 
sion scheme that suggests taking descent steps that can simultaneously reduce both 
data-model misfit terms. Such an approach, however, requires solving a non-smooth 
convex problem at every iteration, which is computationally expensive for a pixel- 
based inversion over the whole domain. Instead, we employ a parametric level set 
(PaLS) technique that substantially reduces the number of underlying parameters, 
making the inversion computationally tractable. The performance of the technique 
is examined and discussed through the reconstruction of source zone architectures 
that are representative of dense non-aqueous phase liquid (DNAPL) contaminant re- 
lease in a statistically homogenous sandy aquifer. In these examples, the geometric 
configuration of the DNAPL mass is considered along with additional information 
about its spatial variability within the contaminated zone, such as the identifica- 
tion of low and high saturation regions. Comparison of the reconstructions with the 
true DNAPL architectures highlights the superior performance of the model-based 
technique and joint inversion scheme. 
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1 Introduction 



In recent years, there has been increasing interest, especiahy with respect to subsurface 
sensing apphcations, in the development of inversion methods that process data from 
highly heterogeneous sets of sensors to obtain a unified characterization of a region of 



space [36,38,40,43,49,68]. These joint inversion techniques are motivated by the idea 
that a variety of modalities, each sensitive to a different set of constitutive parameters, 
combined with appropriate regularization to mathematically relate one set of parameters 
to the others, can significantly improve characterization relative to what can be achieved 
through the processing of individual modalities. The potential of such a joint inversion 
approach has been clearly demonstrated for various modality combinations, e.g., hydro- 
logical and seismic [12], electromagnetic and elastic [3|, electromagnetic and seismic [4l] , 
gravity and seismic 36j, gravity and magnetic [32], as well as magneto-telluric, gravity 
and seismic data (53 . 

In the present work, we consider a joint inversion approach based on a geomet- 
ric parameterization of the problem and employing a new multi-objective optimization 
scheme to combine the data from disparate sensor types. While the method is gen- 
eral, we are specifically concerned with exploring its performance in the context of an 
environmental remediation problem, namely the characterization of a subsurface con- 
taminant source zone based on both geophysical and hydrological data. This work 
is motivated by the problem of remediating sites contaminated by dense non-aqueous 
phase liquids (DNAPLs), such as chlorinated solvents including trichloroethylene (TCE) 
or tetrachloroethylene (PCE), which are used in dry cleaning, degreasing operations, and 
gas production. Successful treatment and management of source zones contaminated by 
such compounds typically is predicated on knowledge of the mass distribution of the 
DNAPLs in the subsurface. 

Partitioning tracer tests (PTTs) are among the most developed methods for DNAPL 
source zone characterization [26|[44] . In this technique, boreholes are used to inject and 
extract tracer fluids. Based on the physics of the problem and a knowledge of observed 
tracer concentration data in the extraction well (and sometimes in boreholes along the 
flow path), an inverse problem is solved to characterize the source zone. There are a 
number of drawbacks with PTTs. For example, this technique only provides spatially 
averaged estimates of DNAPL saturation over flow paths defined by the typically sparse 
distribution of injection and pumping or observation wells |54|. Its reliance on the use of 
boreholes not only may fail to provide full coverage of the affected zone but also increases 
the risk of mobilizing the DNAPL mass, if observation points within the contaminated 
zone are employed. Partitioning tracer tomography is a more recent technique which 
attempts to address the shortcoming of PTT by considering the full transport model and 
characterizing the architecture through a rigorous stochastic inversion. However, this 
method still suffers from the coverage and mobilization risks of PTT. 



In addition to these hydrological methods, geophysical modahties including seismic 
methods IgFI, ground penetrating radar (GPR) and electrical impedance/resistance 



tomography (EIT/ERT) [13,18,34] have also received recent attention as potential tools 
for DNAPL characterization. Less invasive than hydrologic modalities, the utility of 
these methods for the DNAPL problem arises from the contrast in electrical or seismic 
properties of the contaminant relative to the nominal subsurface |7|. Even with these 
methods, data tend to be collected sparsely and one must still solve a challenging inverse 
problem to develop an image of the subsurface. As is typically the case, the need for 
regularization tends to result in images with rather coarse spatial resolution |27|. 

Joint inversion techniques may hold promise in overcoming the challenges discussed 
in the preceding paragraphs. Indeed, a substantial literature has developed exploring 
joint inversion techniques for a variety of related subsurface characterization problems. 
For example, in |47 two-dimensional ERT is performed over various ground transects to 
provide estimates of the structure and water content of the subsurface. These images 
were then interpreted with respect to hydrological investigations of the same basin using 
both tracer methods and groundwater level observations. Researchers have also investi- 
gated coupled inversions based on electrical resistance tomography data and hydrological 
models (e.g., see [38||59|). In these techniques, inversion is performed simultaneously by 
combining the data sets and forming a joint model, instead of constraining the hydrolog- 
ical interpretations via the electrical tomography results or vice versa. A joint inversion 
over these types of data sets significantly improves the estimation and reduces the un- 
certainty encountered with uncoupled techniques. Similarly, GPR data have been used 
along with hydrological data to perform a coupled inversion to identify soil structure |30 . 
In the latter study, the uniqueness and stability of the inversion process was also analyzed 
numerically. 

Although significant strides have been made, the majority of the joint inversion tech- 
niques presented in the literature for subsurface characterization have been for limited 
scenarios (e.g. 2D problems) or simplified physical models ||T2j[30|[47|[58]. Thus, joint 
inversion techniques that can address large scale, complex, fully three-dimensional phys- 
ical domains are highly desirable. Inspired by the idea in |31|, in this work, we propose 
a multi-objective optimization formulation for a fully three dimensional joint inversion 
problem addressing DNAPL source zone characterization. 

The DNAPL problem is one example of a much broader class of inverse problems 
whose primary objective is the identification of a region of interest (source zone, tumor, 
crack in a material sample, etc.) embedded in an inhomogeneous background. For such 
problems shape-based methods, such as active contours and level sets, have received 
considerable attention |24 . Level sets are employed in various applications including 



geophysical and hydrologic inverse problems 12,25,50 . Despite their broad applicabil- 



ity, our prior work in this area has demonstrated that in the case of ill-posed problems 
(such as the one of concern in this paper), traditional level set methods require delicate 



regularization to perform well [5] . The method presented herein employs a new parame- 
terization for joint inverse methods, which is an extension of the alternative parametric 
level set (PaLS) approach, originally developed in the context of single-modality inverse 
methods |5 . In this work we also extend this approach to employ a pair of PaLS func- 
tions to allow for the identification of both high and low saturation regions in the source 
zone, known as pools and ganglia. 

The DNAPL source zone characterization problem addressed here encompasses the 
hydro-geophysical inversion problem illustrated in Figure [T] The hydrologic measure- 
ments for this problem correspond to the observations of contaminant concentration 
collected on a plane (transect) orthogonal to the nominal direction of groundwater flow 
and located down gradient of the source zone. Electrical resistance tomography (ERT) 
data are also collected orthogonal to the flow direction but across planes that intersect 
the source zone. The specific form of this problem is motivated by a number of factors. 
First, the down-gradient transect concentration data are readily available in practice 
and unlike PTT data, their acquisition does not risk any additional mobilization of the 
contaminant. In terms of the geophysical modality, our choice of ERT was driven by the 
fact that there is some contrast between water and DNAPL resistance as well as the ease 
of modeling this modality relative to GPR or seismic data. While GPR may provide a 
more robust signature than ERT for DNAPL in the field ^6^, a rigorous, physics based 
inversion of radar data using the full Maxwells equations is a daunting task. ERT on the 
other hand requires only the solution of Poisson's equation, which is computationally far 
more tractable. 

A preliminary result associated with the problem posed in this paper is reported 
in our recent review |52j; however the material here represents a substantial advance 
over that work. In |52] a PaLS- based approach was used to characterize a simple 
source zone structure whose saturation profile was quite close to uniform. Moreover 
the average saturation value was considered known a priori so that only the geometry 
of the source zone was to be reconstructed. We also assumed perfect knowledge of 
the hydraulic permeability field. As a result, a straightforward scalarized optimization 
formulation performed adequately. In lifting the simplifying assumptions cited above to 
handle more realistic scenarios, we found that a scalarization approach was no longer 
applicable, thereby leading to the development of the multi-objective ideas presented 
here. Additionally, the use of multiple level set functions to identify both pool and ganglia 
is entirely new to this paper. In short, the work in this paper presents a fundamentally 
new variational approach to joint inversion. It extends the PaLS ideas, and demonstrates 
the overall performance of these concepts on a far more challenging form of the DNAPL 
source zone characterization problem than was considered in |52|. 

The remainder of this paper is organized as follows. Section 2 provides a general 
formulation of the problem, discussing the electrical and hydrologic models and their 
relationship to source zone physical properties. Section 3 provides a brief overview of 




Figure 1: DNAPL source zone and associated plume instrumented for hydro- geophysical 
assessment using down gradient concentration observations and cross-gradient electrical 
resistance tomography 

inversion techniques and presents the PaLS representation and the proposed multi- ob- 
jective minimization method. Section 4 discusses implementation of the method and 
presents and discusses some illustrative simulation results. Finally, Section 5 provides 
some general conclusions and implications. 

2 General Problem Formulation 

As we discussed in Section [T} to address source zone characterization, we consider a 
combination of hydrologic and geophysical measurements as illustrated in Figure [T] The 
ultimate goal of this characterization effort is to extract useful information about the 
source zone geometry and the underlying saturation map. In the sequel we provide a 
brief description of each modality and ultimately link them to the proposed joint inversion 
and shape-based characterization techniques. 



2.1 Multi-Phase Transport and Dissolution Model 



Simulation of DNAPL infiltration and subsequent mass dissolution in the saturated zone 
requires the solution of both phase and component mass balance equations. The phase 
mass balance equations are of the form 

^(Pa^5c,)-v-(pa^^(vPa -Pag) =X!E^c.a^ (1) 

which incorporates a modified form of Darcy's equation to represent fluid flow fl]. Here 
pa is the intrinsic mass density of the a phase, (f is the matrix porosity, Sa is the 
saturation, k is the intrinsic permeability tensor of the medium, k^a is the relative 
permeability, /Xq, is the fluid phase dynamic viscosity. Pa is the phase pressure, g is the 
gravity vector, and Eaa'. accounts for the loss or gain of phase mass due to the interphase 
mass transfer of component i from the phase to the a phase. Constitutive equations 
that relate saturation, relative permeability, and capillary pressure, as well as interphase 
mass transfer expressions, are required to close the system The above equation 
neglects intra-phase compositional transformations, e.g., chemical reactions that take 
place within each phase, and since the focus of this study is on the DNAPL mass, 
sorption to the sohd phase is also neglected. 

Within each bulk fluid phase, the spatial-temporal distribution of a component i is 
described by a component mass balance equation written in terms of the mass concen- 
tration of component i in the a-phase (Cf ): 

^1^ (^-CD + . SaiC^^r^ - • VCf) = ^ Eaa^^ . (2) 

Here is the three-dimensional hydrodynamic dispersion tensor for component i in 
phase a p|. The vector v*^ is the a-phase pore velocity computed using a modified 
form of Darcy's law |1 . In this work, the interphase mass exchange of component i 
from the a^-phase to the a-phase (i.e., Eaa') is represented using a linear driving force 
expression |69 . 

For the applications presented here, only two fluid phases are modeled, the aqueous 
phase [a = w) and the DNAPL (a = n), which requires that Sn = 1 - Sy^- Representa- 
tive DNAPL source zone saturation distributions were developed by solving a coupled 
system of equations of the form ([T]) for DNAPL infiltration and redistribution. For the 
plume transport simulations used in the inversion, the redistributed DNAPL is assumed 
immobile and composed of a single component (i.e., i indexing is not needed). Thus, 
solution of only a single (aqueous) flow equation is required and this is coupled to a single 
transport equation of the form Q. Given that sorption has been neglected, dissolution 
is the only interphase mass transfer process considered. 



2.2 Electrical Resistance Tomography 



The ERT model is based on introduction of electrical current into a medium and measur- 
ing the electrical potential at the periphery of the medium to analyze how the electrical 
conductivity is distributed throughout the medium. The underlying partial differential 
equation which relates the potential, i/(x), to the conductivity cr(x) and the electric 
current distribution j(x) is 

V • (aVu) = j , (3) 

with the boundary conditions 

n • Vu = 0, X G , 

n-Vu^Cu = 0, X G Fmix • (4) 

InQ, Tn is a no current boundary corresponding to the air-soil interface modeled using 
a Neumann boundary condition. Over Tmix an infinite half space is approximated using 
a mixed boundary condition by appropriately choosing the function (see [23|[58| ). The 
introduction of current is usually represented as point source dipoles of the form 

i(x) = Jo(5(x-x+)-5(x-x-)) , (5) 

where 6{.) is the Dirac delta function, Jq is a DC current and x"^ are the current electrode 
coordinates. To obtain a denser data set, multiple experiments are performed with 
different electrodes acting as current sources. 

2.3 Petro-physical Relationship 

Petro-physical relationships link cr(x), the conductivity distribution measure of the do- 
main, and in this case, the saturation of the DNAPL, 5^(x) (i.e., a = P(sn))- The most 
widely used petro-physical model is Archies law [8j, which for this two phase aqueous- 
DNAPL system takes the form 

= aa^^^{l-sny. (6) 

Here is electrical conductivity of the aqueous phase, (f is the porosity of the medium, 
a is a fitting parameter, m is a fitting parameter that commonly referred to as the 
cementation index, and q is the saturation index. At large saturations of an electrically 
conductive aqueous phase, the value of the saturation index can be approximated as 
2.0 [28 . Archies Law assumes that the solid and DNAPL do not contribute to the 
electrical conductivity. |42| provided theoretical justification for the form of this model 
using continuum percolation theory. In this paper we use Archie's law, although the 
inversion approach developed is not specific to the selected petro-physical model. 



3 Inversion Strategy 



3.1 Pixel Based and Shape Based Methods 

The goal of most inverse problems is to extract information about a physical property in 
space, p = p(x), using data that are linked to p via a physical model. Consider A4 as the 
computational model that maps p to the data vector d. For simplicity we start with a 
single modality and later extend the notion to more than one model. A straightforward 
strategy to obtain an estimate of p is to minimize model-data mismatch in a variational 
sense: 



p'' = argmin ^ ||d - A^(p) ||r , (7) 

where for a vector u and a symmetric positive definite matrix R 

= u^Ru . (8) 

The matrix R usually contains the noise statistics and a pattern for weighting the data. 
To perform the inversion using conventional approaches, p is discretized over a dense 
grid of pixels in the region of interest and the minimization is carried out to find the 
corresponding pixel values. Given the practical limitations in acquiring dense, rich sets 
of data, many problems of this kind are ill-posed and require regularization. Well known 
regularizations typically take the form of added penalties to the inversion cost function 
to control the amplitude and smoothness of the of reconstruction [4[[66|. 

Shape-based methods are another class of techniques to better pose the problem. A 
shape-based approach proceeds by partitioning the domain of interest into a number of 
zones defined by similar property values. The inverse problem then amounts to determin- 
ing the boundaries of each of the zones, along with a (generally low-order) representation 
for the spatial distribution of the property in each zone. This technique specifically suits 
inverse problems where the main objective is the characterization of an inclusion within 
a background domain. 

The most well known shape-based technique is the level set method [56] , in which the 
shape boundaries are represented via the zero level set of a higher dimensional surface. 
Consider the basic binary case in which the domain of interest, is composed of two 
regions Di and D2, where p(x) = pi in Di and j>(x) = p2 in D2. For a shape-based 
representation, one can characterize both zones using a level set function such that 

6(x) < X e L>2 

and accordingly rewrite the property of interest in terms of pi and p2 as 

piK)=PiH{<t>iK))+p2(l-H{ct>iK))) , (10) 




Figure 2: Level set function and topological flexibility. The orange surface represents 
the level set function (j) and the green plane indicates the zero-level plane. The dark gray 
area of the lower plane is the set of points for which (/>(x) < 0. while the lighter gray 
represent the points for which 0(x) > 0. In three dimensions, this latter set will represent 
the region of the source zone occupied by DNAPL. The bottom frame is obtained from 
the top by "evolving" the level set function according to the velocity field indicated by 
the vectors in the top. Under this motion, the connectivity of the underlying zero level 
set is able to naturally change with an unknown number of components. 

where H(.) represents the Heaviside step function. The scalar anomaly coefficients pi 
and p2 may in general be functions of x, representing some low order representation of 
the anomaly texture in each zone (e.g., see ||45j). The binary case discussed above can be 
generalized to multiple regions by using more than one level set function (e.g., see |12|). 

Minimization of Q for a level-set based property model of the form ( [To| ) is performed 
in an evolutionary fashion. Starting with some initial level set function 00, the function 
evolves to attain a state such that its zero level set best describes the true shape. The 
resulting time-discretized Hamilton-Jacobi type of evolution equation takes the form 

(x) = ,/.W(x)-^W(x) II V0(*)(x) II , (11) 

and is initialized as (/)^^)(x) = 00 (x). In this equation t represents the artificial time in 
the evolutionary process and v^^\:s.) is a normal speed function (shown in Figure [2]). At 
every iteration, the speed function is calculated based upon the sensitivity of the cost 
function to the current shape state |24 . 

An attractive feature of the level set technique is its topological fiexibility. By ap- 
plying a proper speed function, the level set function may deform in a way that informs 
topology changes in the next iteration (note how the orange topology has evolved from 
the top panel to the bottom panel in Figure [2]). Through such evolution, shapes can 
easily merge and split (as shown by the grey shapes in the black plane in Figure [2]) with 
no prior assumptions about the number of connected components in each zone. 



Despite the positive features, there are usually additional complexities associated with 
implementation of level sets, especially for inverse problems |24|. To guarantee proper 
convergence, the level set function should retain a certain form (usually a signed distance 
function). Also, since the acquired speed function only applies to the zero level set of 0, 
speed extension methods to globally deform the level set function need to be applied [55] . 
Moreover, from a numerical perspective, a level set function is still represented in terms of 
discrete grid values (pixels) which potentially increases the dimensionality of the problem. 
Such dimensionality can again cause problems in tackling ill-posed inverse problems and 



pose the challenge of applying traditional or geometric regularizations 10,24|. 

To remedy the problems associated with applying level set techniques to ill-posed 
inverse problems and yet take advantage of their topological flexibility, a recent work 
by Aghasi et al [5j, proposes using a parametric level set (PaLS) function. In the 
PaLS technique, the level set function is parameterized in terms of a parameter vector 
/Xp = [/xi, //2r**5 Mm], with M much less than the number of pixels involved in the dis- 
crete representation of the problem. Its deformations are controlled by changing the 
elements of /ip. This type of parameterization thereby induces a form of regularization 
by parameterization (e.g., see [14,29,45]). In \5\ the authors also present a pseudo-logical 



approach, which approximately models applying set operations on simple geometries to 
form more complex structures. More specifically, considering ^(.) to be a compactly 
supported radial basis function (simply called a bump) they propose a PaLS form as 

M 

0(x,/Xp) = -c+X!«^^A,x^W ' (12) 

where the constant c is a positive scalar close to zero, '0/3.^^. (x) = '0(||/3^(x - Xi) ||) and 
Mp - f^iiXi}ifi is the parameter vector controlling weights, radii, and centers of the 
bumps. Note that c is required to obtain nontrivial results as the radial basis functions 
are themselves exactly zero after a certain radius. As illustrated in Figure [3] this model 
exhibits a pseudo-logical behavior. For example, the sum of two positive bumps having 
comparable size (shown in the top panel of Figure [s] as the orange surface) approximates 
the union operation on their zero level sets (grey shapes in the bottom black plane of 
the figure). Similarly, summation of a positive and a relatively large negative bump can 
approximate the set exclusion operator (Figure [sj bottom). Using this concept, the basic 



algebraic summation in (12) can imply set operations on the support of the bumps and 
make the level set function capable of expressing a large class of geometries. 

By using (/>(x, /Xp) in (|10|) and employing a smooth version of the Heaviside function 
(as discussed more fully inJE])^ the large-scale minimization problem in Q reduces to a 



minimization over the PaLS parameters and the texture parameters, pi and p2 in (10) 




Figure 3: Illustration of the pseudo-logical behavior of the parametric level set functions. 
Top: set union; Bottom: set exclusion 



as 



1 2 

{Mp,Pi,J>2} = argmin - ||d - A^(p(x, /Xp,pi,p2)) || • (13) 



{/Xp,pi,p2}2 



The number of unknown parameters in (13) tends to be much smaller than those 
in pixel based and conventional level set methods. The low dimensionality of the PaLS 
approach makes the inverse problem less ill-posed without sacrificing much in terms of 
flexibility in shape representation. Moreover, it establishes a foundation that supports 
using quadratic minimization techniques such as the Newton methods. These methods 
are faster than gradient descent techniques and robust against the scaling of different 
variables appearing in the minimization [5,33,60|. 



3.2 Joint Inversion and Multi-Objective Minimization 

In Section 2 we presented two different modalities; the hydrological and electrical models. 
In the hydrological inversion we are interested in reconstructing the DNAPL saturation 
values based on the measurements of the contaminant concentration in a down gradi- 
ent transect. On the other hand, in the ERT inversion we seek to extract the electric 
conductivity of the domain given a limited number of potential measurements. Using a 
petro-physical model, the problem can be expressed entirely in terms of saturation and 
inversion can be cast as the solution to the multi-objective optimization problem 

\\\d^-M^{sn)\\l^^ 
st, = arermin < ^ , (14) 

^ ^||ci.--M.(P(..))||^^ 



where indices H and correspond to the hydrological and electrical modalities, and P(.) 
represents the petro-physical relationship that links saturation and conductivity (here, 



Archies law). Using a PaLS shape-based approach as in (13), 5^(x) can be parameterized 
as 5^(x, /x), where /x is a vector containing the PaLS and texture parameters. In this 
case the inverse problem amounts to finding a Pareto solution ^19j for the multi-objective 
minimization problem 

where 

1 2 



GM = ^||d^-A^^(P(5n(x,/x)))||^^ . (17) 



One of the main solution strategies in multi-objective optimization is the scalarization 
approach, for which a single objective function is obtained as a linear combination of the 
underlying costs. For our case, a scalarized version of the multi-objective problem would 
be 

where R^?/ and are fixed at the beginning of the minimization to roughly make the 
two cost terms comparable. In addition to the general problem of balancing the terms, 
scalarization is not always an efficient approach for multi-objective problems (e.g., see 



Section 7 in [31j). In employing gradient-based techniques to iteratively minimize (18), 
there is no guarantee that a descent direction based on Qt simultaneously reduces Qr^ 
and Q^. Moreover, the balance between the two terms may change substantially, as 
the iterations progress, thereby, yielding a solution that in a sense does not maximally 
exploit all of the information in the various data sources. 

For the particular problem of interest here such behavior is in fact observed. Specif- 
ically, DNAPL saturation values may increase without subsequent change to the down 
gradient contaminant concentration, as concentrations along a given flow path approach 
the aqueous solubility. This phenomenon is better demonstrated in Figure Qa). Specif- 
ically in this demonstration, for a fixed volume of DNAPL (here a cube of uniform 
saturation 5), we measure c^, the quasi-steady state downstream concentration of a sin- 
gle point for different values of s. We also have an electric potential measurement, Up^ 
corresponding to a fixed source of current. 

Typical plots of Cp and \up\ in terms of s are shown in Figures IJb) andlJc). Tech- 
nically speaking, once the aqueous solubility of water is reached, for a wide range of 



(a) 




Figure 4: (a) An experiment setup to observe the contrast between the sensitivities of 
electrical and hydrological modalities to the source zone saturation (b) The typical down- 
stream measurement for a single point as the bulk saturation increases (c) The typical 
change in the magnitude of potential measurements as the bulk saturation increases 



saturation values, Cp stays almost constant]^ close to a threshold Cp^^. For the electrical 
measurement, according to the Archie's law and the Poisson's equation, increasing s 
causes an overall reduction in the electrical conductivity of the domain which causes a 
monotonic increase in \up\ as depicted in Figure Qc) . 

This simple experiment reveals that the two physical models would not in general 
maintain a balanced sensitivity to the saturation values in a course of reconstruction. In 

^ After this range, further increase in s would result Cp to start stepping below Cp^^. This reduction 
is due to the fact that as s increases, the permeability to the water phase decreases. At high values of 
s, the water does not flow to any appreciable extent through the contaminated zone and will start to 
flow around it, which decreases the rate of dissolution. As this second phenomenon only happens at very 
large saturations, we only emphasized on the first phenomenon. 



these circumstances the scalarized cost function effectively ignores one data set due to 
larger decreases that can be obtained by "listening to" the other data. While a heuris- 
tic way of approaching this problem is iteratively rebalancing the cost terms, a proper 
convergence may not be guaranteed for such an approach. Furthermore, rebalancing 
the cost terms does not necessarily balance the corresponding sensitivities and still one 
of the costs may be neglected in determining a descent direction. Inspired by the idea 
presented in |31j, in the sequel we present an iterative scheme that the step determined 
in the course of reconstruction would simultaneously reduce both misfit terms. 



3.2.1 A Joint Newton Type Minimization 

Classic Newton techniques provide an iterative procedure to converge to a stationary 
point. More specifically, for a given cost ^(//) initializing the process with fiQ, at the 
k-th iteration a suitable scalar multiple of a descent direction is added to the current 
estimate /X/. to generate the subsequent estimate. To determine the descent direction, a 
quadratic Taylor approximation is considered as 

+ s) ^ ^(/x,) + s^wGifik) + ^s^v^g(fik)s , (i9) 

and to reduce the cost to the maximum extent, Sk is obtained by minimizing g{S) = 

Sk = argmin S^VGif^k) + l^^^^GM^ • (20) 



When the Hessian of the cost, W^Q{iJ,j^), is positive definite, (20) can be solved in closed 
form as Sk = V^G{fJ^k)~^VG{fJ^k)^ which represents the well-known Newton direction. 

In the case of multiple costs Qj{lJi) (where j :'E^9{ here), for every cost we have the 
quadratic approximation 



9j{6) = S'^VGiitik) + \8'^V^QMk)S ■ (21) 

Fliege et al in fSl] discuss that an optimal direction in this case may be acquired by 
minimizing the function q{S) = maxj gj{d). Figure [s] illustrates how a direction acquired 
by minimizing q{S) is the optimal direction to simultaneously reduce both costs Qj{fi). 
Finding this direction requires the solution of the following non-smooth min-max problem 

Sk = argmin max S^yGj{fik) + -^^V^^j(/Xjt)5 . (22) 
S j 2 



Figure 5: Determining a Newton descent direction for multiple costs; the figure shows 
quadratic approximations gi and g2 of the costs Gi and Q2 according to (21), for which 
the corresponding Newton steps are 6g^ and 6g^. Choosing Sq as the ultimate direction 
guarantees that gi{Sq) < and g2{^q) < (i.e., simultaneously reducing both costs Qi 
and Q2). Moreover, as observable in the figure, 5q is in some sense the optimal direction 
that performs this concurrent reduction. 



At the expense of adding an auxiliary variable z, problem (22) can be framed as the 
constrained convex optimization problem 



mm z 



(23) 



Although solving this convex program per iteration is computationally expensive for a 
pixel-based inversion, it suits well into a PaLS framework that /x is of moderate dimen- 
sionality. In other words, the dimensionality reduction that PaLS brings into the problem 



makes taking the Sk steps in (20) tractable and obtain a direction that simultaneously 
reduces both cost terms. 



One of the main limiting assumptions of the multi-objective Newton technique in [31] 
is the positivity of the Hessians, which does not necessarily hold for Qr^ and Q<£. However, 



we can take advantage of the least squares nature of and Gt, in (16) and (17) and 
obtain positive definite approximations to the Hessian. 

Motivated by Gauss-Newton techniques [51], denoting Jj = dMj/d/j, as the model 
Jacobian for j : j^, we have VGj = Jj^Ilj{Aij - dj) and V^Gj 



Jj^HjJj is at least a positive semi-definite approximation to V^^j which by adding a 



Jj Rje/j. 



The matrix 



small positive multiple of the identity matrix, (i.e., Jj^HjJj + Ajl), becomes a strictly 
positive definite matrix. In the Levenberg-Marquardt (LM) algorithm, at every iteration 
an optimal value of A is determined based on trust region or damped techniques [22] . 
Accordingly, the LM algorithm of |51 can be generalized to the multi-objective case, 
where at every iteration by choosing suitable values of Aj, reasonable approximations to 



V Gj matrices are obtained and replaced in (23) to find a descent direction. The details 
of this generalization are provided in [A| 

4 Simulation and Discussion 

In this section the performance of the joint inversion scheme is examined for reconstruc- 
tions of realistic source zone architectures. We first discuss how the models are treated in 
generating synthetic data and then technical details associated with the inversion scheme 
are provided. 

4.1 Data Generation and Underlying Modeling 

Each realization used in this work is representative of a suite of those generated by nu- 
merically simulating a PCE-DNAPL release in a statistically-homogenous sandy aquifer 
medium. Here the University of Texas Chemical Flooding Simulator UTCHEM 9.0 [21] 
was used to solve a coupled system of equations of the form ([T]) for known release and 
flow conditions. Distributions of hydraulic permeability and related capillary parameters 
were produced using sequential Gaussian simulation (SGS) |48|, based upon the general 
characteristics of an aquifer at the Bachman Road glacial outwash site in Michigan |2 ). 
Simulation parameters and boundary conditions are presented in Table [T} Note that, 
for this work, the average permeability and porosity employed are representative of the 
Bachman site, but a larger permeability variance and shorter correlation length were 
applied in the SGS algorithm to produce more variability in the source zone DNAPL 



saturation distributions. The interested reader is referred to 16 and the references 
therein, for more detail related to the infiltration simulations. 

Figure [6] depicts one of the simulated source zones. Here the green isosurface roughly 
delimits the configuration of the contaminated zone, corresponding to a saturation of 1%. 
The brown (darker) spots or surfaces correspond to high saturation zones which exceed 
the residual saturation of these sands (15%) and are typically classified as pools. Also 
shown in Figure [5] (panel b) is a representative 2-D slice of the source zone taken at x =3.97 
m, illustrating the high degree of heterogeneity in the saturation distribution. Note that 
the spatial domain depicted in the figure encompasses all of the infiltrated PCE mass, 
but is smaller than the actual computational domain used to simulate the distribution, 
which extended an additional 5.48 meters in depth. This smaller domain, with the same 




Figure 6: General setting of the problem; a) the saturation profile where the green 
iso-surface corresponds to low saturation values and darker spots correspond to high 
saturations. The electrical setting, water concentration measurements and vadose zone 
are also shown in the figure, b) A slice of the saturation at x=3.97 m showing the DNAPL 
saturation texture in this plane, c) The log valued permeability field corresponding to 
the release. 

grid spacing, was employed for all subsequent hydrologic model computations in the 
inversion modeling process described below. 

The ERT test configuration for the problem is also shown in Figure [5^. The top 
plane {z = 0) represents the air interface (ground surface) where a Neumann boundary 
condition is applied. At the remaining planes, an absorbing boundary condition Q 



is applied to simulate an infinite half space. To account for the potential influence of 
the unsaturated zone in the electrical conductivity measurements, a thin vadose zone 
was added above the saturated domain and beneath the ground surface. The sensor 
placement and configuration are also shown in the figure. In total 130 sensors are located 
on the periphery of the imaging domain. The ten vertical boreholes each accommodate 8 
equally spaced sensors and surface sensors are placed as linear arrays, each with 10 sensors 
filling the gap between a pair of boreholes. To generate a set of electrical measurements, 
32 simulation experiments are carried out within the electrode array. In each experiment 
two sensors act as the current dipole (connected via the lines in Figure [5^) and the 
remaining sensors measure the corresponding electric potential. Current sources are 
placed in six boreholes. The remaining four boreholes, and the surface sensors (indicated 
with lighter color), serve as potential electrodes in all the experiments. The dipoles are 
chosen in a cross medium configuration to make the data more sensitive to the presence 
of DNAPL . Measurement data are generated by a finite difference solution of Poisson's 
equation, over the domain shown, which is discretized into 36x36x36 grid blocks, in the 
X, y and z directions respectively. 

With respect to the hydrological modality, the groundwater flow is assumed to be in 
+x direction. To create the hydrologic observations associated with a particular source 
zone realization, a modified version of MT3DMS [71j, which accounts for rate-limited dis- 
solution and relative permeability effects |[T7[|57] was used to simulate the quasi-steady 
contaminant (PCE) concentrations in a down gradient plane {x = Xmax transect) per- 
pendicular to the groundwater flow (Figure [6]^a)). These concentration data represent a 
single snapshot in time captured during the period of quasi-steady dissolution behavior 
that is characteristic of most DNAPL sites. Here we are less interested in the late-stage 
behavior of a source approaching exhaustion/clean- up but rather are interested in char- 
acterizing the source zone prior to remediation. The permeability field shown in Figure 
was used to generate the "true" DNAPL source zone realization (Figure |6]^a)) and 
its associated quasi-steady down gradient plume transect concentration "measurements" . 
Relevant simulation parameters are presented in Table [Tj Further detail on the simu- 
lation methods may be found in |15|. Note that, in the hydrological model, aqueous 
phase concentrations measurements are assumed known at all grid points in the x = Xmax 
transect. 

We acknowledge that, for this method to become field-practicable, we must eventually 
address the issue that concentration measurements are likely to be sparse. However, 
in this paper we focus on assessing the general performance of the technique before 
considering the more difficult problem based upon sparse data. 



4.2 Inversion 



This section provides details on the inversion strategy. For the transport simulations em- 
ployed in the inversion process, the saturated zone shown in Figure [oj^a) was discretized 
into 26x26x50 grid blocks in the x, y and z directions. For all inversion simulations 
the permeability field was assumed unknown. Thus, the flow, transport, and electrical 
aspects of the inversion assume a uniform, average permeability and porosity. Within 
the inversion, for a given source zone configuration, the groundwater flow field and cor- 
responding contaminant concentrations in the down gradient transect were generated 
using a version of MT3DMS |71 , modified in this work to support parallel computing. 

Archies law was employed for the petrophysical model assuming a uniform porosity 
within the saturated zone. The Archie parameters and porosity shown in Table [2] are 
based upon measurements conducted in our laboratory using Ottawa sands and PCE- 
DNAPL. These parameters fall within the range established within the literature [8| 
[28j . Ottawa sands were used here because they provide a good physical surrogate for 
the Bachman aquifer material (e.g., |2]|62||63j). The electrical simulations employ a 
vadose zone containing an uniform water saturation. Within this vadose zone the water 
saturation was assumed to be 8%. This saturation is consistent with the field capacity 
of the porous media employed in f35|, and represents a second feature (i.e., in addition 
to DNAPL) within the domain that is electrically resistive. Because the vadose zone is 
assumed to have constant porosity and constant water saturation (8%), we calculated the 
effective conductivity shown in Table [2] using Archies Law with the saturation exponent 
reported by |35 . The bulk phase conductivity in all electrical simulations was assumed 
to be 0.05 Sm~^, which is analogous to an aqueous solution of approximately 250 mg/L 
CaCl2. Larger values of the bulk phase conductivity could be employed, but would 
increase the contrast between the DNAPL and aqueous phase, and decrease the difficulty 
of the problem. Use of the relatively low value of electrical conductivity in the aqueous 
phase represents a more stringent demonstration of the methods described herein. 

One of the advantages of the PaLS technique is that it is mesh-less, i.e., the under- 
lying parameters are independent of the discretization. Thus, by using PaLS in a joint 
framework, every modality can have its own discretization method, and the grid points do 
not need to be co-located. This is certainly not the case with pixel-based methods where 
joint inversion on two different set of voxels generally requires some level of interpolation 
to provide a unified representation of unknowns, a process that can be complex. 

An effort was made to account for some of the model and measurement uncertainty 
by applying some random noise to the electrical conductivity and concentration data, 
prior to application of the PaLS method. Specifically, 0.1% additive Gaussian noise was 
added to the electrical conductivity data and 2% Gaussian noise to the concentration 
measurements. Although the current electrical noise level is relatively small compared 
to the total signal, this corresponds to more than 20% uncertainty in the scattered field. 



By definition, the scattered field is the measurement variation caused by including the 
anomaly (in this case DNAPL) in the system [l^. Thus, all the information about 
the inclusion is buried in the scattered field and the 0.1% noise actually adds consid- 
erable uncertainty to the data upon which the inversion is based. Using GPR as the 
geophysical modality may help here, as the contrast we get in the dielectric properties of 
water and contaminants using GPR is larger than the electrical conductivity contrast of 
DNAPL relative to the groundwater |46|. However, fully inverting the GPR data using 
the Maxwell's equations is computationally expensive. Given that the primary goal of 
this paper is to demonstrate the initial utility of our approach for joint inversion, we feel 
that the adaptation required for the processing of field data based on the GPR is a task 
best left to the future. As a result, the noise level in the ERT data is less than what 
would be expected from currently fielded geophysical instrumentation although it is in 
line with electrical systems employed in medical imaging [37[[64| . 

There are currently no explicit adjoint forms to aid extraction of the sensitivity values 
from the hydrologic model. However, use of PaLS significantly reduces the dimensionality 
of the problem, and permits application of a basic finite difference approximation for 
these sensitivity calculations. The Jacobian and model sensitivities for the inversion 
were acquired using the adjoint field technique for the ERT case as described in (5l[6l[6l] . 

To examine the ability of the method to characterize the extent of the source we 
consider a single level set function for the DNAPL representation. In addition, we also 
consider the case where high saturation features are identified within the DNAPL repre- 
sentation. To identify these high saturation features we employ two level set functions; 
one corresponding to low saturation regions (ganglia) and one related to the high satu- 
ration regions (pools). 

4.2.1 Using a Single PaLS Function to Invert for DNAPL Structure 

Based upon the PaLS idea we model the DNAPL saturation as 

5^(x) = 5ii/e(0(x,/Xp)) + 5o(l-i/e(0(x,/Xp))) , (24) 

where Si is a scalar representing DNAPL texture value within the source zone and Sq a 
similar quantity for the texture outside of the source zone and /Xp is the vector of PaLS 
parameters, a, /3 and %. Since the DNAPL saturation outside of the source zone is 
zero. So vanishes and simplifies the inversion such that only the PaLS parameters and 
Si need to be quantified. The function is a twice differentiable approximation of the 
Heaviside step function and e, a parameter which controls the transient width of H^: and 
is a positive number smaller than c (see [5]). 

For the PaLS function, we choose M =45 bumps (Table [2]). The parameter initial- 
izations are intended to be general, and this selection of M is guided by previous work 



that suggested a range of values low enough to ensure computational tractability and 
high enough to capture details [5] . The initial centers Xi randomly chosen around the 
central parts of the imaging box. Values of /3i are roughly set to make the initial shape 
be of comparable size to the domain being characterized. The weights ai are randomly 
initialized to be 1 such that the method starts with a balanced number of positive and 
negative bumps. The texture value Si is initially set to 1%, and is updated through the 
iterative inversion process. Through some experimentation we have found that Si is most 
effectively updated in a sub-iteration using only the ERT data in a classic Newton step. 
The PaLS parameters o;, /3 and however are updated based on the joint data. The 
generic initialization results in the shape shown in Figure [7|^b). 

The result of applying the multi-objective Newton-type (MONT) method is shown 
in Figure [7|^c). The reconstruction visually compares well with the true source (Figure 
[7|^a)) and contains many similar features (e.g., the multi-lobed shape is well represented 
and areas of DNAPL figuring are well captured). Besides the geometric comparison, 
in the simulations we calculate the relative error in the total reconstructed mass. This 
parameter is closely related to how successfully the reconstructed texture parameter 
captures the inhomogeneity of the shape. Based upon this metric, the joint inversion 
presented in Figure [7|^c) is characterized by 5.5% error in the relative mass. To provide 
a more detailed picture of the reconstruction, in Figure [sj^a) we have shown a transect 
of the saturation distribution and the corresponding reconstruction contour. 

In contrast to MONT, reconstruction based upon the scalarization approach for joint 
inversion results in 8.1% relative mass error and a rather poor reconstruction as shown 
in Figure ^d) and Fi gure 8[b). It is worth noting that the success of the scalarization 
approach in providing an acceptable mass error is thanks to the update of the texture 
parameter, 5^, merely based on the ERT modality. 

The reconstruction shown in Figure is based upon the joint inversion of both 
electrical and hydrologic data. Thus, it is important to understand what each data set 
can provide when independently inverted. Reconstruction based on inverting only the 
ERT data is shown in Figures [7|^e) and[8](c). The recovered mass error in this case is 
1%. Reconstruction based on inverting only the concentration data is shown in Figures 
^f) and[8](d) (76.4% relative mass error). The electrical-only reconstruction limits mass 
errors at the expense of shape matching (note the bulkier shape in Figure [7|^e) and 
the loose contour in Figure [sj^d)). In contrast, inversion based only on hydrologic data 
appears to be poor in both relative mass errors and shape matching. This is mainly due 
to the severe ill-posedness of the hydraulic model originating in the thresholding effect 
of aqueous solubility. 

The iterative cost reductions for both the scalarized and MONT algorithms are shown 
in Figure [sj^e) . Both methods are initialized in identical states. In using a scalarization 
approach, the balance between the costs quickly changes and at every iteration the reduc- 
tion of one cost term might be at the expense of increasing the other. The scalarization 



approach converges in about 45 iterations, but ends in a local minimum that is highly 
affected by the hydrologic model (compare Figures [7|^d) and[7|^f)). In using the MONT 
approach both individual cost terms monotonically decrease at every iteration. With the 
MONT approach, descent directions are determined based on the constrained problem 



(23) and the number of iterations before reaching a steady state is about 75 iterations. 

It is important to note that considering a constant hydraulic conductivity brings 
some level of uncertainty into the inversion model. Besides this, there is an added 
uncertainty associated with the thresholding effect of the hydrological model and the 
inaccuracy in the data. Based on these three sources of uncertainty, full matching of the 
hydrologic model with the data does not necessarily result in successful reconstructions 
(revisit Figure [7]^f)). As observable in Figure |8]^e) , in using the scalarization approach 
the ultimate cost value for the hydrologic model is less than that of MONT, originating 
from the dominance of the hydrologic cost in the inversion process. This dominance 
causes an early trapping in a local minimum. In using MONT, however, thanks to the 
mutual control that cost functions have over each other, only the useful portion of the 
hydrologic data is used and by avoiding a full hydrologic data-model matching, an early 
local minimum is prevented. In this case the ERT cost value shows to be less than that 
of the scalarization approach and a more reliable iterative process is achieved. 

4.3 Using Two PaLS Functions for Characterization of Ganglia and 
Pools 

In the previous example a single level set function was used to reconstruct a shape rep- 
resenting the average DNAPL profile. When characterization efforts seek to distinguish 
areas of high saturation, a second PaLS function can be added to create a texture func- 



tion (rather than a texture parameter as described above). A more general form of (24) 
is 



5^(x) = 5i(x)i7e(0i(x,/Xp i)) . (25) 



The main difference between (|24|) and (|25|) is in the nature of 5^. In (24), Si is only a 



scalar while in (25), 5^(x) is a texture function defined within the source zone. This 
texture function is used to classify the source zone texture into low and high saturation 
values. More specifically we represent the texture as 

5i(x) = 5pi7e(02(x, /Xp^2)) + Sg(^l - i7e(02(x, Mp^s))) , (26) 

where Sp and Sg are scalar values representing average saturation values for the pools 
and ganglia. Plugging ([25]) into ([26]) yields 



5^(x) = {Sp - 5^)i7e(02(x, /Xp^2 ) (^1 (x, /Xp^i)) + 5^i7e(0l(x, /Xp^i)) 



In this representation the level set function (f)i characterizes the regions inside and outside 
the source zone. Within the source zone, the function 02 characterizes the ganglia and the 
pools. Therefore this representation is capable of classifying regions as: zero saturation, 
low saturation or high saturation. Here rather than estimating Sp and Sg, for simplicity we 
use fixed values with Sg =0.01 to define the extent of contamination and Sp =0.15 to then 
extract the regions containing pooled DNAPL. These values are arbitrary, but consistent 
with our prior work in defining ganglia and pools [Tt] . Alternative values can be assigned 
to these two parameters, and the inversion will produce the corresponding iso-surfaces. 
PaLS parameters [/Xp i; Mp,2] obtained through the joint inversion technique. 

As a first example in differentiating ganglia and pooled regions we use the same 
electrical and hydrological data sets as the single level set case considered in the previous 
section. In representing (f)i we again use 45 bumps with similar initialization parameters. 
To represent 02 we use 35 bumps, as detailed in Table [2j Initial centers are chosen 
randomly amongst those associated with 01, and the dilation factor is selected to be 
Pi = 2 to result in smaller initial geometries for the pools. The overall initialization is 
shown in Figure [9]^c) where the green iso-surface corresponds to the initial ganglia and 
the brown surface represents the initial pools. 

Visual comparison of the joint inversion results (Figure [9]^c)) with the original satu- 
ration profile (Figure [9]^a)) suggests that the features corresponding to ganglia and pools 
are well reconstructed. The error in recovering the total mass for this reconstruction was 
7.8%. It is interesting to note that use of the two level sets actually decreases overall 
performance of the joint inversion. Thus, the benefits of delineating the high saturation 
regions come at the expense of some accuracy. A transect of the source zone along with 
the underlying contours are shown in Figure [iTj^a). It can be observed that besides the 
good performance in characterizing the source zone, the pool contours are more or less 
accumulated around the denser pooling sites. The pools in this release are very sporadic 
and distributed all around the source zone. To more clearly illustrate the performance of 
the method in characterizing the pools, a second source zone architecture (Figure [To]^ a)) 
was created by simulating the release of a PCE-DNAPL within a different permeability 
field (not shown). For this example although still sporadic, most pooling sites are close 
to the ground surface. Initial PaLS parameters are the same as in previous example 
(Table [2]) though the initialization (Figure [lOj^b)) employed a new set of random centers. 

The resulting source zone reconstruction is shown in Figure [Toj^c) . Besides the suc- 
cessful reconstruction of the ganglia structure, this example more clearly shows the accu- 
mulation of reconstructed pools close to the surface where the actual pools are located. 
Despite the random initialization of the pool locations, in the course of inversion they are 
gradually pushed towards the actual locations. The error in recovering the total mass for 
this example was 2.5%. As before, a 2D transect of the saturation and the reconstruction 
contours are presented in Figure [iTj^b) , that help showing the appearance of the pools 
contour close to the actual pooling sites. 



5 Conclusion 



This paper basically provides a geometric approach to joint inversion. While a specific 
application was considered herein, the formulation and strategy in solving the inversion 
in a joint form is general. In fact, the joint inversion approach that is presented in 
this work can be adopted by other disciplines to combine various physical modalities. 
This is particularly important as many joint inversion techniques simply suggest using a 
scalarization approach for physically incompatible models. 

Application of the proposed technique to the DNAPL characterization problem pro- 
duces reasonable reconstructions of the source zone structure based on the fusion of 
electrical and hydraulic sensors. Using this technique we are able to reconstruct a de- 
tailed picture of the source zone along with useful information about its texture and the 
high saturation regions. In fact, our results suggest how hydrology information can be 
used in a controlled way to "help" improving the geophysical reconstructions. 

Although application of geophysical techniques to 3D problems can be found in the 
literature, use of a fully 3D, multiphase flow and transport model within an inversion 
routine has not previously been demonstrated, due in large part to the computational 
burden and ill-posedness associated with the use of such a model. In the context of joint 
inversion, the problem becomes even more complex, since appropriately combining the 
modalities and making use of different data sets significantly exacerbates the problem. 
Our ability to employ 3D models is directly linked to the low-order nature of the PaLS 
and still its high flexibility in shape representation. Thanks to this low-order representa- 
tion, applying the MONT technique becomes computationally tractable and a controlled 
inversion over incompatible data sets becomes possible. 

The presented work is in fact a proof of concept to what can be done in the future. 
As stated before, a future work would be exploring other geophysical modalities, such 
as GPR, which are capable of providing a higher contrast at the expense of computa- 
tional complexity. Moreover, to reduce the hydrologic model uncertainties, methods of 
estimating the hydraulic permeability prior to (or in line with) the inversion may be 
considered. Ultimately, the successful performance of the method in the challenging and 
rather realistic examples considered, makes it a reasonable technology to be tested with 
real field data. 

A Joint Minimization Algorithm 

The following algorithm is a generalization of the LM algorithm presented in |51| to the 
multi-objective case. For sake of simplicity we consider only two objectives Qr^ and 
(as the case in our paper), however, the algorithm can be easily generalized to handle 
more than two objectives. 



Algorithm 1 A Multi-Objective LM Algorithm 

Vj{ := 2; v<E := 2; £ := sq; % is a small number and can be the machine e. 

/x:= /Xq; 

:= max(diag( J<;/(/x)^R<7/J:7/(/x))); 
Ac := max(diag( J<£(/x)^R<E J<e(/x))); 
found := false; 
while - found do 

Hg{ := J^{fi)^Il^J^{fi) + A^I; 

:= J<e(ij^)^B.<eJ'e(i~J^) + A^I; 
solve (23) with approximate Hessians and H-e 
to obtain a direction 5; 

:= 2(6^^(/x) - g^ifi + d))/(d^(A^d - VGM))'. 
Pe := 2(6^^(/x) - g^^ti + 5))/(5^(A^5 - V6^^(/x))); 
if Pfj{> and p<E>Q then 
/X := /X + (5; 

A<7/ := A<7/max(|, 1 - {2p^ - 1)^); 
A<E := A<E max(|, 1 - (2/><e - 1)^); 
1'^ := 2; — 2; 
found := (\\S\\ < e); 
else 

if p^ <0 then 

A^ := '?;^A^; := 2v^] 
end if 

if p-E <0 then 

A<E := v^X<E] V'E •= 
end if 
end if 
end while 
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Table 1: Source Zone Scenario Simulation Parameters 



Fluid Properties 







Water 


PCE 


Density pa {g/cvc?)^ 




0.999 


1.625 


Dynamic Viscosity (cP)^ 
Compressibility (Pa~^)^ 




1.121 


0.89 




4.4x10"^° 


0.0 


Aqueous Diffusivity (cm^/s)-^ 






8.6x10"*^ 


Aqueous Solubility (g/L)*! 
Initial Saturation 




1.0 


0.150 
0.0 


Spill Scenario 


Spill Volume (L) 






128 


Spill Duration (d) 






400 


Release Rate (L/m^-d) 






0.32 


Pc-Sa-Kra Model Parameters 


Air Entry Pressure (kPa) 






2.809 


Pore Size Index 






2.0773 


Interfacial Tension 








Air /Water (dyn/cm) 
PCE/ Water (dyn/cm) 
Irreducible Water Saturation 






72.75 
47.8 
0.080 


Max Residual Organic Saturation 






0.151 


Reference Permeability (/Lim^) 






19.7 


Matrix Properties 


Variogram Parameters 




Horizontal 


Vertical 


Nugget 




0.333 


0.333 


Range (m) 




4.66 


0.72 


Integral Scale (m) 




1.55 


0.24 


Mean Hydraulic Conductivity, K (m/d)^ 




16.8 


Anisotropy Ratio, k^/kh 






0.5 


Lognormal Transformed K variance 


(<72log(K))t 




1.5 


Applied Hydraulic Gradient (m/m) 






0.01 


Longitudinal Dispersivity, uim (m)^ 






0.30 


Horizontal Transverse Dispersivity, 


'jjp (m)ll 




0.10 


Vertical Transverse Dispersivity, Up 


(m)' 




0.0075 


Median Grain Size, dso (/xm)^ 






295 


Uniformity Index, Ui ^ 






1.86 


Uniform Porosity, cp ^ 






0.36 


Ax (m) (AT^ = 26) 






0.3048 


Ay (m) {Ny = 26) 
Az (m) {N^ = 128) 






0.3048 
0.0726 



^487^20 , ^ 39 



67 



Table 2: A Summary of Underlying Inversion Parameters 



Saturated Zone: 


Size (m) 


7.93x7.93x3.81 




Porosity, ip 


0.36 




(S/m) 


0.05 


Vadose Zone: 


Width (m) 


0.5 




Electrical Conductivity (S/m) 


2.5x 10-^ 


Archie's Law: 


a 


1 




m 


1.4 




Q 


2.0 


PaLS Initialization/Setting: 


Functions and 


Function 02 




M = 45 


M = 35 




e = 0.1, c = OAl 


e = 0.1, c = 0.11 




Xi ' random 


Xi ' random 




A = 0.6 


A = 2 




ai = ±1 


OCi = ±1 




1 

(e) (f) 



Figure 7: Source zone reconstruction using a single level set function: a) original source 
zone structure; b) initialization; c) joint reconstruction of the source zone using proposed 
technique; d) joint reconstruction using a scalarization approach; e) reconstruction using 
only the ERT data; f) reconstruction using only the hydraulic data 
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Figure 8: A slice of the saturation at x = 3.97 m showing the DNAPL saturation in 
this plane and the corresponding reconstruction contour a) Using MONT; b) Using a 
scalarization approach; c) ERT-only data inversion; d) Hydrologic data inversion; e) 
Iterative performance of the scalarization approach compared to the proposed technique 
(MONT). Individual costs corresponding to every technique are marked 




(a)' " ° (b) ' ° ° (c) 



Figure 9: Source zone reconstruction using two level set functions: a) original source 
zone structure with the iso-surfaces corresponding to 1% and 15% saturation values; b) 
initialization; c) reconstruction using the proposed joint inversion technique 




Figure 10: Another source zone reconstruction using two level set functions: a) original 
source zone structure with the iso-surfaces corresponding to 1% and 15% saturation 
values; b) initialization; c) final reconstruction result 



(a) 



(b) 



Figure 11: a) A slice of the saturation presented in Figure [9]^a) at x = 3.97 m; the gangha 
contour is shown in black and the pools contour is in gray b) A slice of the saturation 
in Figure [Tol^a) at ^ = 3.97 m, again showing the ganglia contour in black and the pools 
contour in gray 
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